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Abstract. For the first time a complete set of the most recent direct data on primary 
cosmic ray spectra is used as input into calculations of muon flux at sea level in wide 
energy range Ef^ = 1 — 3 • 10^ GeV. Computations have been performed with the 
CORSIKA/QGSJET and CORSIKA/VENUS codes. The comparison of the obtained 
muon intensity with the data of muon experiments shows, that measurements of 
primary nuclei spectra conform to sea level muon data only up to several tens of 
GeV and result in essential deficit of muons at higher energies. As it follows from 
our examination, uncertainties in muon flux measurements and in the description of 
nuclear cascades development are not suitable to explain this contradiction, and the 
only remaining factor, leading to this situation, is underestimation of primary light 
nuclei fluxes. We have considered systematic effects, that may distort the results of the 
primary cosmic ray measurements with the application of the emulsion chambers. We 
suggest, that re-examination of these measurements is required with the employment 
of different hadronic interaction models. Also, in our point of view, it is necessary 
to perform estimates of possible influence of the fact, that sizable fraction of events, 
identifled as protons, actually are antiprotons. Study of these cosmic ray component 
begins to attract much attention, but today nothing definite is known for the energies 
> 40 GeV. In any case, to realize whether the mentioned, or some other reasons are 
the sources of disagreement of the data on primaries with the data on muons, the 
indicated effects should be thoroughly analyzed. 
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1. Introduction 

A muon component of cosmic rays plays an important role in many fields of astroparticle 
physics. It provides a basis for verification of our knowledge on primary cosmic 
ray (PGR) spectrum behaviour, high-energy hadronic interactions and for solution of 
neutrino physics problems. In view of this, a question of inter consistency of information, 
gained in the last two decades on primary and muon fiuxes, and on high-energy hadronic 
interactions presents great importance. In well known works of Volkova et al. (1979) [1], 
Dar (1983) [2], Butkevich et al. (1989) [3], Lipari (1993) [4], Honda et al. (1995) [5], 
Agrawal et al. (1996) [6] and Bugaev et al. (1998) [7] this problem was outside 
of the consideration, since the principal aim lied in the estimation of the secondary 
lepton fiuxes. Authors relied on the information on primaries and nuclear interactions, 
that was available at the papers writing time. Certainly, these data were incomplete 
and ambiguous, as the consequence the inputs in the calculations also vary rather 
significantly. Nevertheless, the outputs, i.e. muon fiuxes, satisfactory agree with each 
other (see comparisons in [5-7]) and muon experimental data. The given circumstance 
with necessity shows, that discrepancies in the Tiscd PGR spectra in the large part were 
compensated by different approaches to the treatment of nuclear cascades. This, in 
turn, may possibly relate to the fact, that muon flux in most of these works served only 
for normalization of neutrino fluxes and it was flttcd to the muon experimental data 
via adjustment of not precisely known parameters of hadronic interactions. From these 
short remarks it is clear, that all these calculations lack some standardization, as it was 
recently proposed by Gaisser and Honda [8] in respect to the PGR spectra. 

The current situation provides significantly better capabilities for adequate choice 
of PGR model and for simulation of nuclear cascades, and in this paper we have tried to 
take the maximum advantage of this. But, as our analysis has shown, unfortunately the 
information on all components of calculations is still very uncertain to provide a firm 
ground for accurate derivation of muon fiux at sea level and, today, it is following 'from 
top to bottom'. Direct measurements of primary cosmic ray spectra span up to the 
energy of ~1 PeV for protons and up to few hundred TeV/n for other groups of nuclei. 
Most abundant data are collected for EpcR < 1 TeV/n. Here proton spectrum is studied 
in series of recent satellite and balloon experiments with approximately 20% accuracy. 
Higher energy data of SOKOL [9], MUBEE [10], JAGEE [11] and RUN JOB [12] although 
have relatively large statistic and systematic errors, also satisfactory agree with each 
other. At the same time, the measurements of the helium spectrum are much less 
concordant and still differ by almost a factor of 2 for all energies, and this uncertainty is 
especially crucial for evaluation of and u/u ratios. As it is widely known, protons 

and helium nuclei contribute ~ 90% to the nucleon flux at the top of the atmosphere, 
which is relevant for derivation of the muon spectrum at sea level, thus the particular 
emphasis should be given to the accurate description of these spectra. 

To avoid simplification, for the simulation of nuclear cascades we have applied 
widely approved and thoroughly tested Monte-Garlo code GORSIKA [13], that allows 
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to treat hadronic interactions with the use of any of the up-to-date interaction models: 
QGSJET [14], VENUS [15], HDPM [13], NeXuS [16] or DPMJET [17]. Overall 
uncertainty, brought in computations of muon flux by the use of these models, is not very 
significant and considerably decreased in the last decade. So, the difference in p — air 
inelastic cross-sections 'between the models have shrunk from 80 mb to today 20 mb 
in the region of few PeV [18]. Average numbers of muons in cascades from primary 
protons, obtained with the use of these models, differ not more than by 20% [19] for all 
energies of interest. This means, that integral and differential muon fluxes at sea level 
will also differ by a close value, and this discrepancy is the smallest, in comparison with 
the uncertainties of PGR chemical composition and energy spectra for -EpcR ^ 1 TeV/n, 
and the sea level muon data. The latter, for < 100 GeV, are both numerous and, 
on the whole, quite ambiguous (see extensive summary in [20]), but the most recent 
ones of BESS [21] and CAPRIGE [22] are very precise and closely agree with each 
other. They give a basis for checking information on primary light nuclei spectra up 
to the energies ~ 10^ GeV/n. The interval of higher muon energies 10^ — 10^ GeV, 
where the spread of the data on muon intensity is still small and amounts to some 
20%, allows to examine with nearly the same accuracy behaviour of primary proton and 
hehum spectra for the energies, extending to the upper bound of balloon and sateUite 
measurements ~ 1 PeV/n. Evidently, that beyond the PGR data, it is not feasible to 
perform reliable calculations of the muon spectrum. What is more, the inverse problem, 
i.e. reconstruction of PGR ffuxcs from the data on muon spectrum, can not be also 
solved, because for Efj_ > 10 TcV there arc only indirect data, in which muon ffux at 
sea level is derived from results of underground measurements. These data are quite 
contradictory and have large systematic errors, mainly caused by incomplete information 
on rock properties and vagueness in question of prompt muon generation mechanisms 
(see, for example, [7,23]). As a consequence of all this, now there is no possibility to 
make conclusions neither on the preferability of any model of charm generation, nor on 
behaviour of PGR spectra for EpcK > 1 PeV/n. 

As it is stated above, this paper is devoted to the investigation of compatibility 
of the present sea level muon ffux measurements in energy interval E^ = 1 — 10^ GeV 
with the experimental data on PGR spectra for corresponding primary energies 10 — 
10® GeV/n. In section 2 we brieffy review the present data on primaries of H, He, GNO, 
Ne-Si and Fe groups and discuss PGR models, applied in this paper and in the papers 
of the other authors. Basic characteristics of computations are presented in section 3. 
Section 4 is devoted to consideration of consistency of our vertical muon flux, derived 
from the data of direct PGR measurements, with sea level muon experiments in the 
energy range Ey^ — l — 10^ GeV. Since in this section we have revealed a sizable shortage 
of muons with energies, corresponding to that of primaries, studied with the emulsion 
chambers, in section 5 we examine possible effects, that may lead to a distortion of the 
information, obtained in the space and balloon experiments, employing this technique. 
Our conclusions are given in section 6. 
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2. Primciry cosmic ray spectra 

The extensive compilation of the modern data of space and balloon measurements 
plotted vs kinetic energy T, GeV/n is shown in figures 1-3. The most of the data 
are gathered for relatively low energies < 1 TeV/n. But even here some 15 years ago 
the spread in experimental results was at least 100% for all groups of nuclei. In spite of 
series of experiments, performed since then, significant improvement was achieved only 
in the study of the proton spectrum. So, for energies from ~10 GeV, where effects of 
geomagnetic cutoff and solar modulations become negligible, to some hundred GeV, data 
of recent experiments LEAP [24], CAPRICE [25], IMAX [26], BESS [27] and AMS [28] 
are consistent within 20% (see figure 1). Better mutually consistent, within 5-10%, data 
of AMS and BESS, now are considered as the most precise. From ~ 200 GeV to ~ 2 TeV 
there is a gap between magnet spectrometer and emulsion chamber experimental data, 
with the only ionization calorimeter measurements of Ryan et al. [29]. Since at the 
energies, overlapping with AMS and BESS, the given data overestimate proton intensity 
(but not for hehum), probably they should be lowered by 25%, as proposed in [30]. In 
any case, more accurate information on behaviour of the proton spectrum in this energy 
interval is required for reliable evaluation of the muon and neutrino fluxes for the energies 
below 300 GcV. From approximately 2 TcV to 100 TeV SOKOL, MUBEE, JACEE and 
RUNJOB measurements provide rather consistent information on proton flux, though 
with the nearly 20% uncertainty. For understanding of the proton spectrum behaviour 
beyond the 'knee', valuable information is provided by the pioneer measurements of 
the Tibet hybrid experiment [31], carried out at the mountain altitude. These first 
results point at the steepening of the proton spectrum for Ep > 200 TeV and this is an 
essential guideline for the extrapolation of the direct measurements. Spread in the data 
for all heavier nuclei practically for all energies of interest is much larger than that for 
hydrogen, and equals to 100% and more. The difference between the helium spectra 
of SOKOL, JACEE and RUNJOB is of major importance for accurate evaluation of 
secondary lepton ffuxes. 

From the aforesaid it becomes evident, that today no unique flt of all data on 
primary spectra can be given and considerable arbitrariness in the choice of them still 
remains. In such events one may use, for example, upper and lower estimates of all 
primary fluxes, or the flux, giving some average of the experimental data. The principal 
model of the PGR spectrum, employed in our calculations, may be regarded just as 
the latter case. It was obtained in frameworks of anomalous diffusion of cosmic rays in 
fractal interstellar medium [32-35], but here it is considered irrespectively to its validity 
for description of cosmic ray propagation. Altogether in our computations we have 
accounted for 5 groups of primaries: H, He, CNO, Ne-Si and Fe. Mass composition of 
the applied model is fltted to the direct and EAS experimental data on the all-particle 
primary spectrum (flgure 4). As it it seen from flgures 1-3 our model as well satisfactory 
corresponds to the available experimental information on elemental spectra and for high 
energies really gives some average of the data on nuclei, except that for hydrogen it 





Figure 1. Primary proton differential spectrum. Experimental data: [36] Chicago, 
[29] Ryan et al., [37] MASS, [9] SOKOL, [26] IMAX, [25] CAPRICE, [12] RUN JOB, 
[27] BESS, [11] JACEE and JACEE fit, [24] LEAP, [10] MUBEE, [28] AMS, 
[31] Tibet (HD) and Tibet (PD). Dashed line is the spectrum, proposed by Gaisser 
and Honda [8]. Solid line is the spectrum, used in this paper. 



presents rather an upper estimate. For calculation of muon flux it also matters, that 
the data on the proton component become insufficiently statistically provided from 
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Figure 2. Primary helium differential spectrum. Experimental data: [38] Ichimura et 
al. Two lines for the Gaisser and Honda spectrum correspond to their 'low' and 'high' 
helium fits. Other designations are the same as in figure 1. 



~ 10^ GcV, and that sizeable uncertainty gives ±50% spread of the data on hehum 
of JACEE and RUNJOB around the vahics, adopted in our model So, though on the 
whole the situation looks quite well-defined for calculations of vertical muon intensity 
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Figure 3. Primary nuclei differential spectra. Experimental data: [39] HEAO- 

3, [40] Lczniak and Webber, [41] Juliusson, [42] Simon et al., [43] CRN, [44] Minagawa. 
Other designations are the same as in figure 1. 
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Figure 4. All-particle primary spectrum. Experimental data: [9] SOKOL, 
[12] RUNJOB, [11] JACEE, [45] MSU, [46] CASA-MIA, [47] AKENO, [48] Proton, 
[49] Tibet. YAKUTSK, AGASA, Fly's Eye and Haverah Park are cited according 
to [50]. 'Previous measurements' are taken from [51]. Solid line is the spectrum, used 
in this paper. 



up to the energies of several TeV, for the greater energies the existing PGR data are 
very ambiguous. In order to estimate the influence of these data spread on the resulting 
muon flux, we have additionally performed calculations with the use of another two PGR 
models, shown in figures 1-3: first of them presents the JACEE fit of their hydrogen 
and helium measurements [11], and the second is the model, suggested by Gaisser and 
Honda [8]. 

To compute sea level muon fiux we have converted primary nuclei spectra to a 
spectrum in nucleons in standard manner: 

ME) = JpiE) + J2 ME ■ A) ■ A\ 

A 

here Ja{E) is the differential energy spectrum of the nuclei with atomic number A, 
A = A, 14, 28, 56. In figure 5 it is presented, along with our primary nucleon flux (solid 
line), fluxes, applied in calculations of muon intensity in works [1,4-7]. These PGR 
spectra vary both in shapes and values and here it is useful to discuss all of them in 
order to realize, how they correspond to the modern experimental data. First of all, one 
can see that a sizable discrepancy between our, 'Gaisser and Honda' and 'JAGEE fit + 
10%' (here we add 10% to account a contribution of heavier, than H and He, nuclei) 
prevails, and this is a good illustration of the difficulties in the attaining of the unique 
description of the current experimental data. Behaviour of our model requires some 
more comments. Enhancement by ~ 20% of intensity of the proton component in the 
energy range 10^ — 10^ GeV is dictated by the fact, that our first test calculations of 
muon flux with the 'average', close to JAGEE flt, proton flux resulted in a very large 
shortage of muons at sea level (see section 4). In order to smooth this contradiction 
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Figure 5. Primary nucleon differential spectra. [1] Volkova et al.; [53] Nikolsky 
primary spectrum, applied for calculations of muon flux in paper of Bugaev et al. [7] ; 
[4] Lipari; [5] Honda ct al.; [6] Agrawal et al.; [11] JACEE fit+10%. Shaded area is 
the spectrum, proposed by Gaisser and Honda [8], with the lower and upper bounds 
corresponding to the 'high' and 'low' helium fits. Solid line is the spectrum, used in 
this paper. 

we have tested several variants of PGR spectra and eventually adopted the presented 
version. It gives the maximal increase of the nucleon flux, provided the elemental and 
all-particle spectra stay consistent with the experiment. Note, that if the muon data 
did not force us to enlarge the hydrogen flux, then it would reasonably agree with the 
new data of the Tibet AS7 collaboration [31]. Dip of our nucleon and elemental fluxes 
for EpcR ^ 10''"' GeV/n is an intrinsic feature of the anomalous diffusion propagation 
model, and, on the other hand, it corresponds to the numerous indications about the 
change of the PGR spectral index in the 'knee' region. 

Prom the earlier works, the most appropriate PGR models were employed by 
Lipari [4] and Agrawal et al. [6]. The nucleon spectrum from the latter paper almost 
coincides with our one, deviating only from the energy of about 10^ GeV. Mainly, this 
discrepancy is caused by the use of excessive, in comparison with today's, flux of protons, 
obtained from the resuhs of the first 8 JACEE flights [52]. In 1993 JACEE reported, 
that proton spectrum consists of two parts with different slopes and a spectral break at 
Ep = 40 TeV. In paper [6] this break was not accounted for and extrapolation of the 
data for Ep < 40 TeV was applied for the higher energies. In work of Lipari a simple 
power primary nucleon spectrum was picked, which is convenient to use in analytic 
calculations, and no collation with the experiment was presented. Nevertheless, this is 
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a rather appropriate model, closely agreeing with the 'Gaisser and Honda' fit. 

The three remaining spectra, applied by Volkova ct al. [1], Honda et al. [5] and 
Bugaev et al. [7] are quite excessive, compared to what is known today. In the oldest 
of the discussed work of Volkova et al. [1], it is applied very large all-nucleon spectrum 
Jn{E) = 1.9- 10*^£^~^'®^ (GeV-m^ -sec-sr)"^, obtained from EAS data, and no elemental 
composition had been considered. Spectrum, used by Bugaev et al. [7], is taken from 
the work of Nikolsky [53], where it was derived from the analysis of fluctuations of 
muon and electron numbers in EAS. For Ep > 10^ GeV proton intensity in this model is 
essentially overestimated (as a consequence, nucleon spectrum too), in comparison with 
the experimental data, and precisely such behaviour provides a good agreement of the 
resulting muon flux in [7] with the most of experimental data on sea level muon intensity. 
In the paper of Honda et al. [5] the proton flux for all energies is also rather high for the 
following reasons. In the low energy region the given spectrum relies on the compilation 
of Webber and Lezniak [54], i.e. mostly on the Ryan et al. [29] data. But they are 
referred in [54] incorrectly and do not match the original values, overstating them by 
~ 20%. Extrapolation of this spectrum to the higher energies reasonably agree with the 
JACEE 1—8 flights results, but, as it was already said, to some extent overvalues the 
present data. Prom the analysis of these three last models a simple deduction follows: 
their renormalization down to the experimental values should cause a sizable deficit of 
calculated sea level muons. We shall return to this statement again in section 4. 

3. Basic characteristics of calculations 

Simulation of cascade processes in the atmosphere has been performed with the use of 
CORSIKA (v6.00, v6.018). In most of calculations as the model describing hadron- 
nucleus interaction for energies E'lab > 80 GeV QGSJET model is apphed. This model 
is in a good agreement with accelerator and EAS experimental data (see e.g. [19, 55]) 
and provides comparatively high calculation speed. Hadronic interactions with energies 
^lab < 80 GeV were simulated with GHEISHA [56]. 

In order to save machine time and accounting, that contribution of all nuclei, 
heavier than He to muon spectrum at sea level docs not exceed 7-10%, we have used a 
superposition model, considering nucleus with energy as A nucleons with energies 
Ea/A. Validity of such approach for calculation of number of secondaries in EAS at 
sea level is widely known (see, e.g., [57]), however, for verification, we have made our 
own computations of muon numbers in cascades, initiated by primary protons, carbon 
and iron nuclei, for energies 10^'^ and 10^^ eV/n with the QGSJET model. Results are 
presented in figure 6. As it is seen, the agreement between numbers of muons in cascades, 
simulated with realistic nuclei fragmentation and those, obtained in the frameworks of 
superposition models is excellent. This, as well, allows not to make a distinction between 
showers, initiated by protons and neutrons, and further by nucleons we mean protons. 
In figure 7 the relative contribution of all nuclei groups, according to our PGR mass 
composition, to the integral muon fiux at sea level is shown. Since the contribution of 



On inconsistency of experimental data on primary nuclei 



11 



10 : 

10 

1 : 

10-^ - 
10-' - 



Proton X 1 2 ^ 

Proton X 56 v 

Carbon o 

Iron □ 



^primary=10'' GeV/nucleon 



10 lO"" lO'^ 

Eth. GeV 



10'' 
10' 
10 

1 

10-^ 
10-' 



Proton X 12 & 

Proton X 56 v 

Carbon o 

Iron □ 



Eprimary=10^ GeV/nucleon ^ 



10" 



10' 10= 
E,h. GeV 
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He nuclei amounts to ~ 20% throughout the entire energy range, thus mentioned in 
the previous section ±50% spread of experimental results on helium flux around values, 
taken in our model, brings approximately ±10% uncertainty to muon flux at sea level. 
From this point of view, the uncertainty due to ambiguity of the experimental situation 
for heavier nuclei is insignificant. 

To obtain easily the muon spectrum for any model of the primary spectrum, it is 
appropriate to perform calculations of integral muon flux in the following way: 

Emax 

/M(>^th)= J N^{En,> Eti,)MEN)dEN. (1) 
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Table 1. Average number of muons with energy above threshold in shower from 
primary proton. 

Primary energy Ep ~ 10'' GeV 



Threshold energy i?th 


1 GeV 


10 GeV 


10^ GeV 


10^ GeV 




This paper 


1011 


318 


21.0 


0.605 


QGSJET 


[19] 


1085 


316 


20.5 


0.588 




[58] 


1008 


310 


20.9 


0.696 


VENUS 


This paper 


1079 


347 


23.5 


0.679 




[19] 


1150 


349 


24.0 


0.604 




Primary energy 


Ep = 10^ 


GeV 




Threshold energy E^h 


1 GeV 


10 GeV 


102 GcV 


103 GeV 




This paper 


8307 


2292 


132.3 


3.612 


QGSJET 


[19] 


8298 


2207 


124.9 


3.240 




[58] 


8059 


2257 


129.8 


3.332 


VENUS 


This paper 


9629 


2590 


153.5 


3.932 




[19] 


9706 


2653 


155.0 


4.156 



Here N^^En, > E'th) is an average number of muons with energy > Et^ in shower from 
primary nucleon with energy Ej^, Jj^{En) — differential primary spectrum converted 
to spectrum in nucleons, E'max — energy, which provides the calculation accuracy 
Ifj,{> E'th) ~ 0, 1%. Most of results, presented below, obtained at £'max/-E'th = 3 ■ 10^. 

For each threshold energy average muon numbers were computed for 20-25 different 
primary energies (see figures 8, 9) with accuracy, generally better than 5% . Vertical 
lines in these figures show the areas of primary energies, giving 10%, 50% and 95% 
contributions to the integral muon flux in case of power primary spectrum Jn{E) = 
1.9 ■ lO'^-E'^^'^^ (GeV • m^ • sec ■ sr)~^. It is seen, that primaries with energies within 
i^th — SOOE'th on 95% determine muon intensity at sea level. Totally, with the use of the 
QGSJET it was simulated about 3 • 10^ showers. Besides, for examination of sensitivity 
of muon spectrum to hadron-nucleus interaction model, series of calculations with the 
VENUS model have been carried out. 

For the verification of our results, we have performed a comparison of numbers 
of muons in showers from primary protons, obtained by us, with those from papers of 
CORSIKA authors [19,58] for two primary energies Ep = 10^ 10^ GeV (table 1). On 
the whole, we see good agreement between the data, presented in this table. Some 
difference in muon numbers for i^th = 10'^ GeV may be explained by a statistical error, 
since simulation of 500 showers, made in the above mentioned papers for each primary 
energy, provides approximately 8% and 4% mean square deviation for energies Ep — 10^ 
and 10^ GeV correspondingly. 
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Figure 8. Number of muons 7V^(£p, > Eth) with energy above threshold, generated in 

shower from primary proton with energy Ep. Vertical lines show the areas of primary 
energies, giving 10%, 50% and 95% contributions to the integral muon flux in case of 
power primary spectrum Jm{E) = 1.9 • 10^E~'^-^^ (GeV ■ m^ ■ sec • sr)~^. 



4. Results and discussion 

Obtained integral and differential muon spectra for our model of primary spectrum are 
presented in figures 10, 11. It is appropriate to split the consideration of the situation 
into three energy intervals. 

Efj^ G [1 — 10^] GeV. In this region our calculation of differential muon spectrum is 
in a good agreement with the most recent and accurate data, obtained by BESS 1995,97- 
99 [21] and CAPRICE 1994, 1997 [22]. This, along with the availability of unambiguous 
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Figure 9. Relative contribution of primary energy region from E'th to Ep to the 
muon integral flux at sea level Iij,{> E^h) versus Ep/E^]^. Vertical lines have the same 
meaning as in figure 8. 



data on primary spectrum and on behaviour of hadronic cross sections, once again 
evidences in favour of correctness of the appUed calculation procedure. Selection of 
BESS and CAPRICE experiments is motivated by the fact, that measured there muon 
spectra have small statistic and systematic errors and excellently correlate with each 
other. In these experiments measurements were made with superconducting magnet 
spectrometers. Earlier data, obtained with iron magnet spectrometers differ from each 
other not only in values, but also in shapes of measured spectra, that may be due to 
the influence of some improperly accounted systematic errors. More detailed discussion 
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of this question may be found, for example, in [20,21,59,60]. 

Approaching to 100 GeV, the calculated muon flux becomes deficient, conforming 
only to the CAPRICE data. Two circumstances in this connection are worth noting. 
The first is that already from these energies information on muon intensity is not definite 
enough. The second is that muons with i?^ > 100 GeV are most effectively produced 
by the interactions of primaries with EpcR > 1 TeV/n, measured in space and balloon 
emulsion chamber experiments. To this moment we shall pay more attention in section 5. 

G [10^ — 10^] GeV. Experimental data on differential proton spectrum 
for corresponding primary energies 10^ - 10^ GeV of SOKOL, MUBEE, JACEE and 
RUNJOB show good interconsistency, but from 10^ GeV they become less definite for 
technical and natural (low fiux) reasons. Experimental data on muon component have 
smaller errors and agree within ~ 20%. In this region along with direct measurements 
on the surface (MARS [61], Nottingham [60,62], L3 [63,64]) there are results obtained at 
underground installations, in which the sea level muon spectrum is reconstructed from 
the 'depth-intensity' curve and in other ways. The latter are the data of the Baksan 
neutrino observatory (BNO) [65,66], Artyomovsk scintillation detector [67], MSU [68], 
Prejus [69], LVD [70], MACRO [71], KGF [72,73] and CosmoALEPH [74]. Most of them 
are grouped within fit hmits, given by MACRO, and comprise a well correlated data 
set. Comparison of the obtained differential and integral muon spectra with these data 
does not support reasonable expectation 'PGR data fit^nuclear cascades^muon data 
fit'. The deviation of our results from the data of 'MACRO zone' is about 30-40%. This 
fact requires revising all steps of our calculation in order to realize how real the problem 
is and what may be done to overcome it. First of all, one should analyze an adequacy of 
nucleon cascades description. As it is shown in the previous section, our computations do 
not contain methodical errors, and total statistical + interpolation + integration error 
is less, than 5%. Prom figure 11 it is also seen, that the maximal difference < 10% in 
predictions of the QGSJET and VENUS hadronic interaction models is also practically 
out of significance in this situation. Since the VENUS model was initially chosen by us 
as the model, providing one of the largest muon numbers in shower (see [19]), hence 10% 
should be regarded as the utmost possible increase of muon intensity for any interaction 
model in comparison with the QGSJET. 

In this connection the results of the earlier works, not indicating any problems with 
the description of the muon data, deserve some closer consideration. We analyze here five 
works [1,4-7]. The collation of the PGR spectra with the experimental data, performed 
in section 2, allows to split these papers into 2 groups: papers with excessive [1,5,7] 
and appropriate [4, 6] fiuxes. In order to reveal differences in hadronic interactions and 
calculation technique, we have made computations of muon fluxes for the same PGR 
spectra, as used in the discussed works. Results are presented in figure 12. From this 
comparison it becomes evident, why no muon deficit was discovered earlier. In addition 
to overvalued nucleon fiuxes, in the papers from the first group the muon yield is larger 
or, at least, comparable with that of QGSJET and VENUS. Muon production in the 
papers with appropriate primary fiuxes is 20% — 30% higher, than obtained by us. 
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Figure 10. Differential muon spectrum at sea level. Experimen- 
tal data: [61] MARS 1975, [62] Nottingham 1968, [60] Rastin 1984, 
[75] MASS 1989, [22] CAPRICE 1994,97, [59] OKAYAMA 2001, [76] Bate- 
man 1971, [21] BESS 1995, 97, 99, [63] L3 1993, [64] L3 2003, [74] Cos- 
moALEPH, [70] LVD 1998, [66] Baksan 1992, [67] ASD 1985, [68] MSU 1994, [69] Fre- 
jus 1994, [71] MACRO fit 1995. QGSJET and VENUS arc the present work 
calculations with the corresponding interaction models. Prompt muons spec- 
tra: [77] — RQPM (recombination quark-parton model), QGSM (quark-gluon string 
model); [78] — VFGS (Volkova, Fulgione, Galeotti and Saavedra model). 
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Figure 11. Integral muon spectrum at sea level. Experimental data: [60] Rastin 1984, 
[62] Nottingham 1968, [75] MASS 1989, [61] MARS 1975, [72, 73] KGF 1990, 

[65,66] Baksan 1990, 1992. QGSJET and VENUS are the present work calculations 
with the corresponding interaction models. Prompt muon spectra: [77] — RQPM, 
QGSM; [78] — VFGS. 
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and this difference is sufficient to mask the muon shortage. Besides, PGR flux from the 
work [6] becomes excessive from ~ 10 TeV, as wc have already pointed in section 2. The 
data in respect of muon production of the latter two papers closely agree, apparently 
because in both of these works essentially the same authors participated. In [6] also 
an explanation of relatively large muon yield may be found. In this calculations, for 
simulation of hadronic interactions TARGET code had been applied, but with the 
enhancement of pion and kaon production. This resulted in increase of contribution of 
kaon decays to muon flux with E^> 1 TeV up to a very large value of 50%. Of course, 
uncertainties in description of hadronic interaction prevail, but, as our calculations show, 
discrepancies in muon yield at sea level between the widely recognized and extensively 
tested models, included in GORSIKA, are at the level of ±10%. These models are 
evidently free from many drawbacks, peculiar to the models from the discussed works, 
and more thoroughly developed. Hence, likely the largest part of the discrepancy should 
be attributed to underestimation of primary nucleon flux, but not to the incorrectness 
of the simulation of the nuclear cascades in the atmosphere. 

There is no need to repeat arguments from section 2 in favour of adequacy our PGR 
model, since it may be regarded just as a valid representation of the current experimental 
situation. We should only repeat, that our proton spectrum does not present a fit of 
the experimental data for Ep = 10"^ — 10^ GeV, but even on ~20% overestimates them. 
This was done to smooth the discussed problem. If directly input in calculations fits of 
primary H and He spectra, obtained by JACEE, with addition of 10% contribution of 
heavier nuclei, or 'standard' spectrum, proposed by Gaisser and Honda, then resulting 
integral muon flux will be even lower, than ours (see table 2). 

There is a possibility to compensate some of the muon deficit by increase of intensity 
of heavier nuclei in our model. But to get 10% addition to muon fiux it is necessary 
to rise hehum fiux 50% or to double aggregate fiux of GNO, Ne-Si and Fe groups. 
This mechanism is rather ineffective, because limitations, set by the data on all-particle 
spectrum thus would require to lower fiuxcs of other, lighter species. Summarizing, we 
see, that today for this energy region the situation is determined rather rigidly and 
the stated problem of muon deficit can not be resolved without some extraordinary 
assumptions. Necessity in rise of all PGR spectra as high, as possible to diminish the 
muon shortage, indicates, that measurements with emulsion chambers systematically 
understate fiuxes of the primaries. A large spread in the data on PGR nuclei with 
A > A provides some more ground for such speculations and demonstrates necessity in 
further improvement of experimental technique. This question is closely related to the 
appropriate description of the characteristics of nuclear interactions and in more details 
is considered in section 5. 

Efj^ G [lO'^ — 3 • 10^] GeV. Analysis of this interval in the present situation looks 
vain for the following reasons: 

i. Obviously, muon deficit takes place for these energies too, and its nature and 
quantity are unclear; 

ii. There is no experimental data on chemical composition of PGR spectrum for 
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Figure 12. Ratio of integral (a) and differential (b) muon fliixes, obtained with 
CORSIKA, to the muon fluxes calculated in papers [1] (Volkova et al.), [7] (Bugaev et 
al.), [6] Agrawal et al., [5] Honda et al., [4] Lipari. The CORSIKA results are for the 
same PGR spectra as used in these papers. 



Table 2. Ratio of the integral muon intensity, calculated with the use of the present 
paper PGR model, to the intensities, obtained with the JAGEE PGR fit [11] and the 
model, proposed by Gaisser and Honda [8]. 
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energies higher than 10^ GeV, that is why, as in our case, one has to use only theoretical 
guidelines, describing the origin and transport of cosmic rays in interstellar medium. 

iii. Behavior of hadron-nucleus cross-sections for E\^\^ > 10^ GeV, despite the 
decrease of divergence between predictions of different models, remains also not enough 
studied, that results in additional uncertainty in calculations of the muon spectrum. 

iv. All muon component measurements for the given energies are obtained from 
results of underground experiments and contain quite large systematic errors, that causes 
their poor mutual agreement. Dominant sources of the mentioned errors are such factors 
as incomplete information on chemical composition of overburden rock and necessity to 
input in calculation a ratio X of prompt muons to tt, i^-muons as a function of energy, 
which is not only unknown, but should be determined from the experimental data. 

Taking aforesaid into account, it seems impossible to make any definite conclusions 
on behaviour of PGR spectrum and on prompt muon contribution on the basis of sea 
level data in this energy region. Though from the comparison with the experimental data 
the VFGS (Volkova, Fulgione, Galcotti and Saavcdra) model [78] of charm generation 
looks as the most preferable, however, here it is needless to stress all unreliability of this 
deduction. 
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5. Possible reeisons of systematic PGR flux underestimation 

In this section we intend to discuss some problems, that are inherent to the study 
of high-energy cosmic rays spectra with the application of the emulsion chambers 
(EC). In addition to all factors, hsted below, one can easily add some more, that 
experimenters encounter during processing of the EC data, including, for example, not 
negligible experimental bias. We are mostly concentrated at the factors, depending on 
the description of nuclear interactions, and in this aspect the situation may eventually 
turn out similar to that with the calculations of sea level muons. The necessity of 
further improvements in the experimental technique is clearly recognized in connection 
with the serious disagreements between the current data on different groups of nuclei, 
and short history of space and balloon EC experiments gives many examples of how the 
reported data changed with the gain of experience and statistics. Further investigations 
will inevitably provide information, that improves our knowledge, and, possibly, will 
differ from it. 

The most detailed description of the experimental technique is presented by the 
RUNJOB collaboration [12] and we shall mostly rely on their data, bearing in mind 
that the procedure is generally the same in all EC experiments (MUBEE, JACEE 
and RUNJOB). The key moment lies in the determination of the energy, transferred 
to the electromagnetic component ^^e cascades, initiated by the interaction 

of a primary particle within EC. This energy is related to the initial energy via 
partial inelasticity coefficient k (we shall denote so k^ to avoid further confusion with 
the spectral index 7). The spectrum of primaries is obtained from the spectrum of 
electromagnetic cascades (EMC) with a simple shift in the energy scale by the value 



here f{k) is a distribution function of k and 7 is the spectral index of primary spectrum 
J{E) ~ E~^'^~^^\ The given statement, at the least dating back to as early as 1962 [79], 
holds true, provided f[k) does not depend on Eq in a wide energy range. Usually 
independence of C{k,'~f) on energy is argued from the evidence, that total and partial 
< k > inelasticities are constant. But this may prove to be incorrect (see, e.g [80]), 
since < k'^ > essentially depends on the contribution of the large values of k to the 
distribution f{k), and even though < k > is constant, 'tail' of f{k) may change with 
Eq. Besides, large body of information on behaviour of the total inelasticity, presented 
in the literature, points at its slow growth with energy. Another important moment, is 
the difference in < A; > between predictions of various interaction models. In comparison 
with FRITIOF, applied by RUNJOB, VENUS gives ~ 10% lower value of partial 
inelasticitiy < /c > in a single p-A, Ylc-A collisions [12] (other experiments do not provide 
information neither on the code used, nor on the single interaction characteristics). 
Note, that 10% correction to the energy conversion from ^ E^ to Eq for power PCR 
flux J{E) ~ E'"^-^ is equivalent to ~ 30% correction to the intensity. Thus, at first 
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sight, the use of the VENUS provides not only the largest number of muons at sea level, 
but, as well may lead to an increase of the measured proton and helium fluxes. In fact, 
the estimation of the influence of the discussed effect on the results of EC measurements 
is much more complicated. Firstly, one needs to evaluate total energy, transferred to 
EMC in successive interactions of secondary particles, i.e. from the whole shower, not 
from single interaction, and to obtain distribution function f{k) for this fraction k of 
the initial energy. Secondly, this distribution for k is rather broad, so the conversion 
to £'0 for individual shower is impossible and it is necessary to derive effective value 
of C(A;, 7), allowing to get PGR spectrum from the EMC one. So, simple knowledge 
of differences in < /c > in the first interaction does not allow to make straightforward 
qualitative estimate of changes in the final energy shift value. The analysis becomes 
even more complicated, if to account, that the choice of the interaction model strongly 
infiuences the previous steps of the experimental data processing. Namely, it affects 
the determination of ^ and detection efficiency. For evaluation of both of these 
values a complete simulation of nuclear cascades in EC is required. For example, in the 
RUNJOB experiment the 'actual' EMC energy ^ E^^i^ue is obtained from the estimated 
with the 7-ray core method one E^^^sti via direct Monte-Carlo simulation of showers 
with the FRITIOF code. In the JACEE experiment [11,81] determination of cascade 
energies also involves complete Monte-Carlo calculations of the transition curves, used 
for calibration of the ^ E^^ derived from the direct electron counting or from the X-ray 
film densitometry. Analogous calculations are required for evaluation of the detection 
efficiency (see detailed description in [12]). 

Thus, deviations between the interaction models in hadron-nucleus cross-sections, 
energy spectra, multiplicities and phase space distributions (this is important for 
estimation of the most energetic in the lab system and back-scattered in CMS numbers 
of gamma-quanta. VENUS predicts the largest quantity of the latter, compared to 
the other models [19]) of seconday particles, in fluctuations of the energy, transferred 
to EMC, etc., may signiflcantly influence the interpretation of the EC measurements 
results. The mentioned arguments give enough ground for performing of thorough 
analysis of sensitivity of the EC data to the use of various interaction models, which 
is still missing. This as well would allow to make fully consistent estimations of the 
cosmic ray fluxes, i.e. with the application of the single model, both at the top of the 
atmosphere and at the sea level. 

Besides, it is possible that some data distortion may come from changes or anomalies 
in characteristics of hadronic interactions at very high-energies. At present, these 
questions arc widely discussed in connection with the data of EAS experiments, also 
employing EC technique. 

Another hypothetical reason, which may cause systematic errors in processing of 
the EC experiments data, could be the presence of sizable fraction of unusual component 
in PCR. This may be regarded equivalent to the 'change' in hadronic interactions and 
all aforesaid about that effect holds true in the given case. In this connection, rising 
interest to the primary antiproton component deserves closer attention. 
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Figure 13. The primary antiproton to proton fluxes experimental data. Direct 
measurements: [82] CAPRICE98, [83] CAPRICE94, [84] BESS95,97, [85] IMAX92, 
[86] MASS91, [87] Bogomolov et al., [88] Buffington et al., [89] Golden et al. Upper 
estimates of the p/p ratio from the ground-based experiments: [90] — solid arrows; 
dashed arrows: [91] Tibet, [92] L3+C, [93] MACRO. 



As it is appropriately supposed, antiprotons present natural compound of PGR, 
produced mostly in interactions of nuclei (dominantly protons) with interstellar medium. 
Other, exotic p generation mechanisms, result in softer spectra and they are out of 
interest for the energies, we are dealing with (for more information, see references in 
caption to figure 13 and e.g. [94]). The experimental study of p component began 
in 1970's from the first balloon measurements (Golden et al. [89] and Bogomolov et 
al. [87]), and today the direct data cover region from few MeV to approximately 40 GeV 
(figure 13). Higher energy region is scarcely studied and is marked by few indirect upper 
estimates of p/p ratio. One of them by Stephens (1985) [90] is derived from the analysis 
of sea level muon charge ratio. Soon we shall present elsewhere such calculations with the 
employment of improved since then knowledge on PGR fluxes, hadronic interactions and 
muon charge ratio. Other three p/p estimates by Tibet [91], L3+G [92] and MAGRO [93] 
collaborations are made with the use of the Moon and the Sun shadow effects. The 
idea of this approach is the following: positively and negatively charged particles are 
deflected by the Earth's magnetic field to the opposite sides, and when the Sun or the 
Moon, during their transit over detector, block particles, deficit of secondary particles 
should be observed, if antimatter presents in PGR along with the matter, at the both 
opposite sides of the real Sun/Moon positions. Since these experiments failed to detect 
'antimatter shadow', they could report only upper limits of p/p ratio to be equal 22% at 
10 TcV [91], 27.5% at 20 TcV [93] and 15% for muon sample with > 65 GeV (primary 
energy ~ 1 TeV) [92]. Theoretical predictions are not able to clarify the situation, since 
there is no generally recognized propagation theory, reproducing satisfactory even low 
energy data on p/p ratio together with all totality of the data on elemental nuclei 
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abundances, secondary to primary nuclei ratios, and fluxes of positrons and gammas. 
Standard difi^usion/reacceleration models fail to do this, requiring introduction of speciflc 
assumptions, generally working for description of some kind of observations, but leading 
to contradictions in other cases. Among such assumptions one may readily call necessity 
in 'artificial breaks in diffusion coefficients and in the primary injection spectrum' [94] 
to match both p/p and secondary to primary nuclei ratios, suggestion about harder 
interstellar proton and/or electron spectra to explain well known 'GeV excess' in diffuse 
gamma ray flux, etc. (for more information see, e.g. [94-97] and references therein). 
The detailed discussion of these problems is beyond the scope of our paper and we shall 
return to it elsewhere. 

In any case, it is clear, that there is no flrm ground to state neither the presence of 
significant fraction of p's for high energies, nor the reverse and that actual p/p ratio is 
still to be determined via direct and various kinds of indirect observations. Thus, one 
can not completely exclude probability, that among events, identified in EC as protons, 
sizable fraction of antiprotons presents. As it is already said above, it is difficult to 
estimate the influence of the differences between p-A and p-A interaction characteristics 
on the results of PGR measurements with EC. It is known that antiprotons deposit more 
energy into 7's, than p in a single interaction [19], but, we should repeat, the complete 
analysis of changes, produced in the measured flux requires to perform evaluation 

of many parameters, that is possible to do only via direct Monte-Carlo simulations with 
respect to the specific experimental conditions. 

We have made similar estimates for the total muon fiux at sea level. As we have 
elucidated, actually all primary nucleons in our calculations are protons. Replacing 
them by antiprotons, we have computed integral muon flux for Eth = 10^, 10^, 10^ GeV 
and have not found any noticeable changes, in comparison with the calculations for 
protons. In other words, total vertical sea level muon flux turned out to be insensitive 
to the fraction of antiprotons in PCR and, possibly, the same may happen with the 
EC results. Of course, this is not true for //+ and //~ spectra, and ratio, their 

calculations may help to determine presence of p's in PCR. 

In conclusion, we would like to present an estimation of the changes in nucleon flux, 
required to match the muon experimental data. In figure 14 we show, along with initially 
applied spectrum, its part for ii^ > 2 TeV 8% shifted toward the higher energies. From 
the magnet spectrometers region to the EC energies we use an interpolation. So, muon 
spectra, given in figures 15, 16, are calculated with the use of the combined nucleon flux, 
shown in flgure 14 with open circles. These spectra still present rather a lower estimate 
of the muon data and may be approximated via the following formula (energy is in GeV, 



with parameters, given in table 3. 

The used nucleon spectrum is just a model spectrum, since in practice, correction 
to the PCR flux hardly may be constant and should depend both on the energy and 



El = 10 GeV) 



Si,(E) = ^£;(B+Cln(S/Ei)+D 



(2) 
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Figure 14. Initially applied nuclcon spectrum (solid line) and its part for the region, 
relevant to the EC measurements, 8% shifted in the energy (dashed line). Circles show 
the spectrum, used to get the lower estimate of the sea level muon flux, presented in 
figures 15, 16. 

Table 3. Parameters of approximation (2) of muon spectra at sea level. 

Integral spectrum (cm^ • sec • sr)-i 

E, GeV A B C D F 

<100 5.8392-10-3 -8.4964-10-1 -2.5331-10-1 6.7652 - 10-^ 3.3739 - IQ-^ 
> 100 6.5584-10-3 -9.1783-10-1 -2.3051-10-1 1.3452 - 10-^ -3.3109-10-^ 
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atomic mass of the projectile. But it is easy to note, that required 8% shift is significantly 
smaller, than the quoted accuracy of the energy determination in the EC experiments 
and, on the other hand, to some part it may be attributed to the effects, discussed in 
this section. 



6. Conclusions 



For the first time a complete set of the most recent experimental data on PGR nuclei 
spectra has been used in calculations of muon flux at sea level in a wide energy range 
from 1 to 3 • 10^ GeV. The adequate description of cascade processes in the atmosphere 
has been provided by the apphcation of widely approved code GORSIKA along with 
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Figure 15. Differential muon spectrum at sea level. QGSJET and VENUS — present 
work calculations with the corresponding interaction models for the primary nucleon 
spectrum, shown in figure 14 with open circles. Other designations are the same as in 
figure 10. 
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Figure 16. Integral muon spectrum at sea level. QGSJET and VENUS — present 
work calculations with the corresponding interaction models for the primary nucleon 
spectrum, shown in figure 14 with open circles. Other designations are the same as in 
figure 10. 
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QGSJET and VENUS models. This makes our results free of the possible errors, caused 
by a simplified approach to description of hadronic interactions. 

The main and a rather unexpected conclusion, which our investigation has led to, 
is that the data on primaries conform to the muon data only for rather low energies 
up to several tens of GeV. From these energies a deficit of muons becomes evident, 
corresponding to a deficit of primary nucleons with £^pcr ^ 1 TeV. In the earlier works 
of the other authors this problem was masked by the uncertainties of the data on PGR 
fiuxes and hadronic interactions, existed those days. In a surprising way, refinement and 
gain of the information have not approached us to the possibility of the more accurate 
muon fiux derivation, but, on the contrary, have revealed that our today's, improved 
notions are in some part erroneous. This, by the way, makes obsolete attempts to 
standardize the lepton fluxes calculations in respect to the primary spectra, since, as 
our consideration shows, in the large part the problem lies in underestimation of PGR 
intensity. In our point of view, to remedy the situation with the least consequences, 
the re-examination of the EG experimental results with the application of different 
interaction codes is required. As we have noted, in this connection VENUS possesses 
with two positive features. In comparison with the other models it predicts a larger 
number of muons at sea level, and in the single hadronic interaction, less fraction of the 
energy is deposited into electromagnetic component, but the influence of the latter effect 
on the measurements results is not so simply to analyze, since it may be compensated 
at the different stages of the EC data processing. Another fact, that must be taken 
into account, is the possible presence in PGR of sizable fraction of antiprotons. If the 
emulsion chambers data would prove to be insensitive to the given effects, only then 
the problem should be attributed to some uncertainties or anomalies in the high-energy 
hadronic interactions. In any case, our calculations have demonstrated, that today the 
fit of the sea level muon data can not be derived directly from the PGR data, and the 
reasons, leading to this puzzle, look unclear. 
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